import matplotlib.pyplot as plt
import numpy as np

d = np.array([15]) * np.ones([18, 1])
f = 50
lmbda = 1500 / f
N = 18
n = np.arange(0, N, 1).reshape([-1, 1])

# 信号
theta = -10
theta = theta / 180 * np.pi
x = np.exp(1j * 2 * np.pi * np.sin(theta) * n * d / lmbda)
Rx = np.mat(1000 * np.dot(x, np.transpose(np.conj(x))) + 1 * np.eye(N))

# # 波束形成
theta = np.arange(-90, 90, 0.1).reshape([1, -1])
theta_0 = 10
theta = theta / 180 * np.pi
theta_0 = theta_0 / 180 * np.pi

a = np.mat(np.exp(1j * 2 * np.pi * np.sin(theta) * n * d / lmbda))

a_theta_0 = np.mat(np.exp(1j * 2 * np.pi * np.sin(theta_0) * n * d / lmbda))
w = Rx.I * a_theta_0 / (a_theta_0.H * Rx.I * a_theta_0)

B = w.H * a
B = np.abs(B) / np.max(np.abs(B))

fig, ax = plt.subplots(1, 2, sharex=True, sharey=False, figsize=(7, 2), dpi=300)
plt.subplots_adjust(wspace=0.35, hspace=0.05)

ax[0].plot(theta.reshape([-1,1])*180/np.pi, 20 * np.log10(B.reshape([-1,1])), '--')

# 方位估计
theta_0 = np.arange(-90, 90, 0.1).reshape([1, -1])
theta_0 = theta_0 / 180 * np.pi

sigma = []
for i in range(theta_0.shape[1]):
    a_theta_0 = np.mat(np.exp(1j * 2 * np.pi * np.sin(theta_0[:, i]) * n * d / lmbda))
    sigma.append((1 / (a_theta_0.H * Rx.I * a_theta_0)))
    
sigma = np.array(sigma).reshape([-1, 1])

ax[1].plot(theta_0.reshape([-1,1])*180/np.pi, 10 * np.log10(sigma))
plt.show()


